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We introduce a simple general method for finding the equilibrium distribution for a class of 
widely used inexact Markov Chain Monte Carlo algorithms. The explicit error due to the non- 
commutivity of the updating operators when numerically integrating Hamilton's equations can 
be derived using the Baker-Campbell-Hausdorff formula. This error is manifest in the conser- 
vation of a "shadow" Hamiltonian that lies close to the desired Hamiltonian. The fixed point 
distribution of inexact Hybrid algorithms may then be derived taking into account that the fixed 
point of the momentum heatbath and that of the molecular dynamics do not coincide exactly. 
We perform this derivation for various inexact algorithms used for lattice QCD calculations. 
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I. INTRODUCTION 

The first algorithms used to simulate lattice QCD with dynamical quarks fell into three classes: Monte Carlo [ll, 0] , 
Langevin [3||i and Molecular Dynamics (microcanonical) i5(] methods. The first two involve taking many essentially 
independent small random steps, leading to a dynamical critical exponent z ~ 2 typical of Brownian motion, where 
this exponent relates the autocorrelation time of the Markov process to the correlation length of the physical system. 
The steps are small either in the sense of making only a small global change in the gauge field or by making a larger 
change locally, such as updating a single link variable. The third suffers from the problem of not being ergodic in 
general, and therefore not necessarily generating the correct distribution of configurations; nevertheless, although the 
updates are also built out of many small steps these steps are correlated so as to lead to a dynamical critical exponent 
z — 1, and it can thus be considered as a "large step" algorithm. The last two algorithms also depend upon a step 
size parameter 5t, and are only "exact" in the limit 5t —f 0. 

It was then realized that by combining Molecular Dynamics with a momentum refreshment heatbath an ergodic 
large step "Hybrid" algorithm resulted, although it still suffered from step size errors in the equilibrium distribution. 
Shortly afterwards it was found that by combining such Hybrid updates with a Metropolis Monte Carlo acceptance 
test these step size errors could be completely eliminated [7|: this is called the Hybrid Monte Carlo algorithm. Both 
the idea of combining Molecular Dynamics with momentum refreshment ^ and that of correcting the Langevin 
algorithm by a Metropolis test [9j had been introduced previously in other fields. 

The introduction of the Hybrid Monte Carlo algorithm did not end the use of inexact algorithms for two reasons: 
firstly some people believe that the volume dependence of the cost of the Hybrid Monte Carlo algorithm outweighs the 
advantage of its exactness, and secondly they wanted to simulate with what are now called "rooted staggered quarks" 
for which there is no explicit local form for the action, but instead the square- or fourth-root of the corresponding 
determinant is required. A way of doing this with step size errors of 0{St) was introduced with the R algorithm 
p^ . Only recently has the RHMC algorithm provided an efficient exact alternative. 

The purpose of the present paper is to analyze the inexact algorithms mentioned above, as well as several other 
interesting variants. The approach we take is new, namely we view all the algorithms as a combination of a Molecular 
Dynamics trajectory with a momentum refreshment; from this point of view the Langevin algorithm is just a single 
step Molecular Dynamics trajectory. Although in principle, any numerical integration scheme could be used for the 
Molecular Dynamics integration in practice all the algorithms use symmetric symplectic integrators (or closely related 
integrators). Such integrators have several remarkable properties, such as being reversible and area preserving, which 
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are required for the exact Hybrid Monte Carlo algorithm. We make use of the remarkable property that all symplectic 
integrators exactly conserve a "shadow Hamiltonian" close to the desired one for small enough step sizes to produce a 
uniform and simple asymptotic expansion for the equilibrium (fixed point) distribution of the corresponding inexact 
Markov chains. 

We then turn to the class of inexact algorithms that use a noisy estimate of the force due to the fcrmionic deter- 
minant. We introduce the pseudofermionic x algorithm as an aid to establishing that the R algorithm, which may be 
considered as an interpolation between the x and Rq algorithms, has only 0{5t^) errors in its equilibrium distribution. 
We also establish that one our original motivations — namely trying to find an variant of the R algorithm with errors 
falling as a higher power of the step size — is almost certainly doomed to failure as the unwanted noise contributions 
to the error are necessarily positive. 

This paper is organised as follows: in f|TT] we use the Baker-Campbell-Hausdorff (BCH) formula to derive the 
form of the "shadow" Hamiltonian that is conserved by various symplectic integrators. A brief description of Hybrid 
stochastic algorithms is given in ijnil and the resulting fixed points of these algorithms are derived in ijlVI We then 
turn our attention to noisy algorithms, specifically as applied to fermion theories, in Jvl For the reader's convenience 
we collect as appendices derivations of some important results that are rather difficult to find in an accessible form 
in the literature. 



II. SYMPLECTIC INTEGRATORS 



We are interested in finding the classical trajectory in phase space of a system described by a Hamiltonian H{q,p) = 
T{p) + S{q) = ip^ + S{q). It has been known for a long time that the leapfrog integration scheme has many desirable 
properties, and higher-order generalizations have been discovered several times in different fields [H, [IE [IE [I3j 
[isl . [iqI . [20| . [21I [2^ . [23! . The basic idea of such a symplectic integrator is to write the time evolution operator as 



exp 



dt 



exp 



exp 



exp T < — 



dp d dq d 
dt dp dt dq 
dH d dH d 
dq dp dp dq 



where H is the Hamiltonian vector field. In the second step we have made use of Hamilton's equations, and thus 
implicitly the fundamental 2-form^ u) = dq A dp. 

From the structure of our Hamiltonian, namely the fact that the kinetic energy T is a function only of p and the 



potential energy S* is a function only of q, it follows that exp ~tS' {q)-^ 
Let us write 



so that by Taylor's theorem 



Q = T'{p) 



d_ 
dq 



3*Q 



and exp 



and 



P EE -S'{q) 



d_ 

dp' 



are trivial^ to evaluate. 



(1) 



f{q + tT'{p),p), 
f{q,p-tS'{q))- 



(2) 



^ For simplicity of presentation we shall consider here only the case where the fundamental symplectic 2-form \s u) = dq /\ dp. For gauge 
theories we need to use ui = — ^Di '^(Pi^i) where the 0i are the left-invariant Maurer— Cartan forms on a Lie group manifold [24j . but 
this generalization is straightforward. 

^ For gauge fields the p lie in a Lie algebra and the q in the corresponding Lie group, so the evaluation is straightforward if not entirely 
trivial. 
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then from the BCH formula (|B3|) we find that the QP symplectic integrator leads to the time evolution given by 

= (exp [(Q + P)5t + i[Q, P]5t'' + ©(.Jr^)] 
= exp[r((0 + P) + i[Q,P]5T + O(5r2))] 

= e^H' =e'(Q+P) +0{5t). 

In addition to conserving energy to 0{5t) such symplectic integrators are manifestly area preserving. 

The BCH formula tells us more than that there is an area-preserving approximate integration scheme for Hamilton's 
equations: it tells us that for each symplectic integrator there exists a "shadow" Hamiltonian H' close to the original 
one, which is exactly conserved. For the QP integrator the actual trajectories through phase space are integral curves 
of the vector field^ 



H' = Q + P+\[Q,P]5t+^[[Q,[Q,P]]-[P,[Q,P]]]5t' ~^^[P,[Q,[Q,P]]]5i 

+^{ -4 [P,[Q,[Q,[Q,P]]]] -6 [[Q,P],[Q,[Q,P]]] +4 [P,[P,[Q,[Q,Pm 

-2 [[Q,P],[P,[Q,P]]] - [Q,[Q,[Q,[Q,P]]]] + [P, [P, [P, [Q, P]]]] ]5t^ + 0{5t^) 

^ -p{S'S"' + S"^)^5t'' +^U[-p^S^^'^ +p{(!>S"^ + iS'S"')) ^ 



Y2" \ s s 



dq 



+ - p'{iOS"'S" + 65(4)5') - 125'(25"' + S'S'")) ^y^^ + ^((Jt^) 

from equation (|B3p . In Appendix El we show that the vector field H' is a Hamiltonian vector field, that is there exists 
a Hamiltonian H' such that 



dp dq dq dp' 

and that H' may be obtained explicitly by replacing the commutators in the BCH expansion (|B3p of ln(e"^e"^) by 
Poisson brackets'* of the form 

{A B} = _ dAdB_ 

dp dq dq dp ' 

so 



H' = H+^{T,S}Sr+^^(^{T,{T,S}}-{S,{T,S}})ST^~j^{S,{T,{T,S}}}dT' 

+ rh{ -4 {5,{T,{r,{T,5}}}} -6 {{r,5},{r,{T,5}}} 
+4 {S,{S,{T,{T,S}}}} -2 {{r,5},{5,{r,5}}} 
- {r,{r,{T,{r,5}}}} + {5,{5,{5,{r,5}}}} )sT' + oi6T') 

= H+ \pS' 5t + ^ (p^S" + S"^) 5t^ + -^pS'S" St^ 

+ (-/^^^^ + P^QS'S'" + 125"') + 125''5") St^ + 0{6t'). 



^ Observe that H' is linear in and because all but that first term in the BCH expansion are commutators. 

Oq op ^ 

* The Poisson brackets define a Lie algebra since they are manifestly antisymmetric and satisfy the Jacobi relation {A, {B, C}} + 
{B, {C,A}} + {C, {A,B}} = 0, q.v., Appendix^] 
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Note that H' cannot be written as the sum of a p-dependent kinetic term and a g-dependent potential term. This 
means that any hope that one could exactly "invert" this relation to find a Hamiltonian whose approximate integral 
exactly conserves H is forlorn. 

A. Symmetric Integrators 

It is immediately apparent that we can do better than this by using a symmetric symplectic integrator (jB4|) which 
just gives the familiar PQP leapfrog scheme for which 



(exp [{P + Q) 5r - ^ ([P, [P, Q]] + 2[Q, [P, Q]]) 5t^ + 0{5t'>) 
\P + Q)-^. {[P, [P, Q]] + 2[Q, [P, Q]]) + 0{5t^) 



r / 5t 



= exp 

^e^(P+Q) +0{5t'). 
The PQP integrator exactly conserves the Hamiltonian 

H' =H+j^ {2^25" - S'^} 5t^ + ^ {-/^W + 6p2(yy" + 25"') - 35''^"} St" + 0((5t^), 

whereas the QPQ integrator conserves 

H' = H+^ {-p'^S" + 2^'} St^ + ^ [ip^S^^^ - 24p2(35'5"' + S"^) + 965''^"} St" + 0((5t<'). 

It is possible to construct symplectic integrators of arbitrarily high order; some aspects of this theory are discussed 
in Appendix ICl 

III. HYBRID STOCHASTIC ALGORITHMS 
A. The Hybrid algorithm 

The Hybrid algorithm [2^, [2^, so called because it is a hybrid of the Langevin Q and Molecular Dynamics 
algorithms, constructs a Markov process on a "fictitious" phase space consisting of the field variables of interest (the 
coordinates) and a set of corresponding "fictitious" momenta. It should be emphasized that these momenta have 
nothing to do with the momenta which occur in the field equations of motion and canonical quantization relations; 
they are just introduced to define a classical dynamics of the fields in a new "fictitious" time dimension. If we are 
considering a four dimensional field theory, for example, then our "fictitious" dynamics takes place in a new fifth time 
dimension. To this end we introduce a Hamiltonian H{q,p) and a fundamental symplectic 2-form as in section|TTl We 
then iterate three Markov steps, each of which has the distribution e^^ as an approximate fixed point and which, 
when taken together, are ergodic.^ 

• The first such step is momentum refreshment whereby the fictitious momenta are chosen from a Gaussian 
heatbath.^ 

• The second step is to integrate Hamilton's equations using an approximate integrator for some length of ficititious 
time r, and then to flip the momenta 



U{t) 



p 



It has the distribution e~^ as an exact fixed point, and it thus approximately preserves the desired Hamiltonian 
H. If the integrator is symplectic then this step is area-preserving, and if it is also symmetric then the step is 
reversible (this is why we incorporate the final momentum flip). 



In some cases we should choose trajectory lengths from some random distribution for this to be true. 

A g eneralization of this is the "partial momentum refreshment" 27] step used in the second-order Langevin [28l.l29j or Kramer's algorithm 

[lallil,!!!. 
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• The third step is a momentum flip, F : p i-^ —p, which obviously preserves the Gaussian distribution of momenta. 
Since this step is immediately followed by a momentum refreshment it can be ignored in all cases but Kramer's 
algorithm [11 [li, [13, [sil, [3| . 



B. Langevin Algorithm 

The Langevin equation [1,[3^ is just a special case of the Hybrid algorithm for a single leapfrog integration step with 
PQP ordering: in other words we use the symplectic integrator e.^'^PI'^^^'^Q e^'^PI'^ using the notation of equation ([T]). 
This operator maps ((7,p) ^ {q',p') where q' = q + pdT — ^S' {q)5T^ , which is the familiar Langevin equation when 
the the Gaussian distributed momenta p are written as Gaussian noise rj and the time step as e = 5t^ : 

q' = q- iS"((?)e + r;Ve- 

The PQP symplectic integrator has the amusing property of giving rise to a simple closed form Langevin equation, 
but there are many other variants too. For example, the QPQ symplectic integrator ^^'^QI'^^^'^p ^^'^Ql'^ leads to the 
equation 



the lowest-order PQP Campostrini wiggle (q.v., Appendix [C|) gives 

^ [(6(3 + ^ + V^)S"'S' - 3(4 + 4^ + 3^)^"') p 
-(3 + ^ + ^)5(4)p3] St^ + 0(^t6) 

and the second-order PQP Campostrini wiggle gives 

q' = q + p5T-^S'Sr^-^S"p5T^ + ^[S"S'-S"'p^]dT^ 



(3) 



1 

' 120 



35"'5' + 5"')p-5(V 



St" 



1 

'720 



-{3S"'S' + S"^)S' + (5S"'S" + eS^^^S'^ p2 _ ^(5)p4j ^ Q(^^7) 
which agrees with the exact evolution operator e^"^ ^ up to 0((5t^). 



IV. EQUILIBRIUM DISTRIBUTIONS 

We now want to address the question as to what fixed point (equilibrium) distribution is produced by the Langevin 
[33I [33 , [35} and Hybrid algorithms. The existence of such a fixed point distribution and its uniqueness follow from 
their ergodicity. We are, of course, looking for a fixed point distribution of q alone, and not of both q and p. 

Since in general the momentum dependence of H' is not exactly Gaussian we do not expect either or 
to be an exact fixed point of a Hybrid Markov process. As we expect the equilibrium distribution to be close to the 
desired one, e~'^, it seems reasonable to parameterize it as e~^'^+^'^). The condition that this is a fixed point is 

g-(s(9')+AS(g')) 

■^^^-(S(,)+AS(,)) r dpe-iP'S{q' -q") 



= J dqdpe-"^'^ P^-^'^^^'^d{q' -q"), 



where ^ ^ ^ ^^^'^ P ) ^^"^ phase space point reached at the end of a trajectory of length r. 
To solve this equation we change variables to {q",p"), whence we obtain 

-(s(.q')+ASiq')) 



J dq" dp" (det C/,)-^e-(^+^^)°^"5(g' - q"). 
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If ?7 is area-preserving then its Jacobian det (7* = det ^^^^^j^ = gtrinC/, _ -|^^ and if it is reversible then (7 ^ = FoUoF 
where F is the momentum-flip operation F : p i—^ —p. 

In the general case we introduce the operation 6 : H. t-^ flo U^^ — floFoUoF, which is zero if U is reversible, so 

{H + AS) o = (H + AS) o F o U o F + 6{H + AS). 

Next we note that H is an even function of p, namely H o F = H , and AS does not depend on p, so AS o F = AS 
also, hence 

{H + AS) oU-^ = {H + AS)oUoF + S{H + AS). 

In terms of the operator 6:i^i-~^ftoUoF — ftoF, which measures the change in some quantity fl over a trajectory,^ 

{H + AS)oU-^ 

= {H + AS)oF + 5{H + AS*) + ~5{H + AS*) 
= H + AS +{d + 5){H + AS). 

Therefore 

g-(s(9')+AS(5')) 

= / dq" dp" e-(^+'^'S')e-('5+^)(^+'^'S')-trinC/.^(-^/ _ ^//-j 

= g-(s(9')+AS(9')) J^y/g-ip""g-(5+J)(H+AS)-trlnC/.^ 

and thus we obtain the condition 



-((5+(5)(_ff+AS)-trln(7, 

Since H is extensive so is 5H, and thus we can show order by order in 5t that AS" is extensive too 



= 1- (4) 

p 



A. Langevin Algorithm 

For a reversible and area-preserving integration scheme, such as a symmetric symplectic integrator, the equilibrium 
distribution must satisfy 

^g-^(H+A5)^ =1. (5) 

For the PQP Langevin integrator we have 

6H = ^ {-p^S'" + SpS'S"} 6t^ + ^ { 

+3p^i2S'S"' + S"^) - 35''5"} 6t^ + ©(^r^), 

and 

SAS ^pAS' St + {^p'^AS" - ^S'AS'} St^ + 0{St^). (6) 
If we expand the integrand of equation ([5]) we obtain to leading non-vanishing order in St 

{SH + 6AS)^ ~ 0, 

and thus 

^{S'AS' - AS"}St^ 

_ 1 + s"S'^ + 5(4) - 2S"'S'} St^ 

+0{5t''). 



An extra momentum flip is included for convenience. 
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This differential equation has an integrating factor of e , 

{[AS' + l{S"S'-S"')dT^] e~^)'~0, 

so 

AS' + 1{S"S' ~S"')6t^ ^ Ke^' 
for some constant K. As we are only determining the asymptotic expansion of AS* we set K — 0, whence we obtain 

AS' ^ i {S'" - S'S"} 5t\ 
Up to an constant, which is fixed by the normalization of e~('5+^S')^ this gives 

AS^l \2S" - ^''} 6t^ + 0{6t'') 

for the equilibrium distribution corresponding to the PQP Langevin equation. 

It is most important to realise that this is only an asymptotic expansion, and the exact solution for AS may 
also involve terms which are exponentially small in the integration step size St. One way to understand this is 
to observe that equation is a Gaussian integral over the momenta p, and that the domain of integration must 
include momenta p ^ 1/6t. In particular this means that the statement made in that the shift in the equilibrium 
distribution corresponds only to the addition of irrelevant operators, and that therefore the St errors can be ignored 
if one computes quantities in the continuum limit, is erroneous because of the existence of these subleading relevant 
contributions. 

If we write the asymptotic expansion as AS — X]n>2 ^"^n then the preceding calculation is easily extended to 
find that the next term satisfies the equation 

~ ± (^S"'S'^ + 2S"S"' - 25(4)5' + 5(5)) . (7) 

While this does not give a closed-form expression for A54 for arbitrary 5, it is obvious that if 5 is a polynomial in 
q then A54 and all the other A5„ are too. In general the fields q have multiple components, and equation ([7]) then 
becomes 

(A54)i ^ ^ (SijkSjSk + 2SjkSijk — 2SijjkSk + Sijjkk) 

with the obvious notation. 

For the QPQ integrator the shift in the equilibrium distribution is 

A52 - i(5''-5"), 

A5; = (245"'5' + 85"'5'' - 205"5"' 
-45(4)5' -5(5)) ; 

whereas for the PQP Campostrini wiggle (q.v. , Appendix [C]) the leading shift satisfies 

A5; - ^ (2(4 + 3^ + 2^)5"'5' 

-(6 + 4^ + 5^)5"'5'' - 6(4 + 3^ + 2^) 5" 5"' 
+2(6 + 5^ + 5^)5(4)5' - (6 + 5^ + 5^)5(5)) . 

B. Hybrid Algorithm 

For the Hybrid algorithm we cannot derive an explicit formula for the shifted equilibrium distribution in general, 
but we can easily see from equation pi Ap that the leapfrog integrator conserves a Hamiltonian H' which differs from 
H by terms of 0{St'^). This means that SH = 0{St^), and thus we deduce from equation ^ that i5A5 ~ 0{St^). 
Unlike the Langevin case considered previously the change in A5 over trajectory has no reason not to be of the same 
size as A5 itself, and thus we find that A5 ~ 0{St'^). 
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V. NOISY ALGORITHMS 



One of the main advantages of the Langevin and Hybrid algorithms is that they can be used for some non-local 
actions where an unbiased stochastic estimate of the force^ S' can be computed relatively cheaply, with (Y}'{q)),f — 
S'{q) where the average is over some "noise" r] which is chosen independently for each step. 

It is helpful to think of the dynamics as that of a system evolving in the presence of some fixed background noise 
field; we shall leave this dependence on rj implicit except where necessary. The noisy force estimator corresponds to 
the discrete mapping e*^ of equation ([2|) becoming 



„tp 



and we thus find that the symmetric QPQ integrator conserves the Hamiltonian H' = H + AH where^ 



AiJ = S - 5 + ^ 



y E" + 21]'' 



This means that the change in H over a trajectory {qi,Pi) i— > {qf,pf) is 

SH = H{qf,pf) - H{q,,p,) = [H{qf,pf) - H{q,,p,)] - [H'{qf,p}) - H'{q,,p,)] 

= [H{qf,pf) - H'{qf,pf)] - [H{q,,p.,) - H'{q,,p.,)] = -AH{qf,pf) + AH{q,,p.,) 



= -SAH = SS 

= iS'~^')pST+^ 



j_ 

24 



T.''' -S'T.' -p'iH" -S") 



Si 



+^ [p3(45"' - 3S"') + 6p(2E"I] - S'E" - 25"!]')] St^ + 0{St^) 

where the last line is for a single leapfrog step where we have used equation ([3]). 

If we consider the quantity e^*^ averaged over the noise 77 for a single leapfrog step we obtain 



(8) 



Observe that the coefficient of (1 — p^)St^ is proportional to the variance of the noisy estimator E', and thus can only 
vanish if the force is computed exactly. Thus, regardless of the order of the integrator used there will always be an 
0((5t^) leading error for a single MD step. 

A. Noisy Langevin Algorithm 

If we now use equation ([5]) to compute the equilibrium distribution for the noisy Langevin algorithm we must 
average equation ([S]) over a Gaussian distribution for p, thus we find that 



where 



A 



-85'(E'E"),, + 4(1]"')^ - 25"(E''),, 



+85''5" + 4(E"'l]')„-5'5"' 



We do not necessarily require a noisy estimator S for the action S itself to implement the equations of motion, but we will require one 
if we wish to construct an "exact" algorithm, q.v. i|V Fl 
® There are two contributions to AH: the first due to the use of an approximate integrator in the presence of the background noise field, 
the second due to the noise itself. 
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Since, for a noisy force equation ([6]) is now of the form 

SAS = pAS' St + ^{p'^AS" - T,'AS')St'^ + 0{St^), 

we immediately find that AS = 0{St^), and more specifically it satisfies the equation e'^ (e~'^AS'')' = —2ASt^. 

B. Noisy Hybrid Algorithm 

In the case of the noisy Hybrid algorithm (for which there are t/5t steps per momentum refreshment) we recall that 
the evolution is to be averaged over independent noise for each integration step. Since the momentum distribution is 
no longer Gaussian after taking a leapfrog step, the leading order term is not cancelled as it would be if we only did 
one MD step per trajectory. The fixed point condition now becomes 

L-S{H+AS)\ 

\ I p,rii,...,riN 

where SH = SHi + • • • + SHn- Thus the leading error of the noisy Hybrid algorithm is AS ~ 0{tSt), as each step 
contributes an error of 0{5t^) by equation ©. A way of reducing this error to 0{tSt'^) is the R algorithm 10], which 
we shall discuss in WEI 

Note that when we construct a Hybrid algorithm with a noisy force estimator there is no reason to expect there to 
be a nearby conserved "shadow" Hamiltonian, as averages of non-linear Poisson brackets will not be correct. Thus if 
the average error per step is 0{5t^) we expect the errors at the end of a trajectory of t/St steps to be 0{tSt) and 
not just 0{St). 



C. Pseudofermions and the $ Algorithm 



We now turn our attention to the specific case of gauge theories in the presence of dynamical fermions. We start with 
the probability distribution for the gauge field U with the quadratic fermion contribution integrated out in favour of the 
fermionic determinant, P{U) cx e"'^'^'-'^-' det A1(C/), where Sq is the pure gauge part of the action and M is the fermion 
kernel. For the purpose of this discussion we ignore the pure gauge contribution to the action since this is a simple 
local quantity whose force can be computed exactly, and shall focus on the awkward fermion determinant. As usual we 
may replace the determinant with an integral over a complex pseudofermion [l| field $ and write the joint probability 
distribution of the gauge and pseudofermion fields as P{U, $) oc exp {- [So{U) + ^■fA^-i^] } = e-^-«(^'*). We have 
taken the fermion kernel to be = JVP M representing two flavours of fermion with Dirac operator M to allow for a 
simple implementation of the pseudofermion heatbath. 



1. $ Algorithm 

The gauge field U corresponds to the variable q used in the general discussion before, and we will evolve it along a 
classical trajectory in the presence of a fixed pseudofermion background field <&. We introduce a conjugate momentum 
field p in order to define a Hamiltonian H{U,p), P{U,p) cx exp {— [ip^ + Scs{U)] } = e~^'^'^^; the action Scs{U) 
takes the role of the potential in the Hamiltonian, and of course depends implicitly on the background pseudofermion 
field. Note that the pseudofermion force term is given by 



The $ algorithm of Gottlieb et al. \lu\ is identical to the Hybrid algorithm described in gHl] with the addition of 
pseudofermion refreshment from a Gaussian heatbath before each MD trajectory. With regards to fermions, the $ 
algorithm is restricted to describing Nf degenerate fiavours of fermions, where iVf is the number of fermions described 
by the operator A4. 

The $ algorithm uses a QPQ integrator, as the evaluation of the pseudofermion force acting on the gauge fields 
required for the P step is only evaluated once in this case. 
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According to the general arguments presented in ijlVI the fixed point distribution of the $ algorithm must satisfy 
equation Performing an asymptotic expansion of this in powers of St, and observing that SH — 0{6t^) for any 
trajectory length r, we deduce that 5/S.S ~ 0(5r^). We thus have that AS* ~ 0{6t^), so the $ algorithm is accurate 
to this order. 



2. X Algorithm 

A slight variation of the <I> algorithm is the x algorithm. The difference is that in the latter the pseudofermion 
heatbath refreshment is performed before every MD step as opposed to before each MD trajectory in the former. We 
might expect the error per step to be 0{5t^), leading to an error per trajectory of 0{t5t). However, the error per 
trajectory is in fact 0{t5t'^),^^ a proof of this will follow from that of the R algorithm to be given in W El 

We note in passing that the pseudofermion force for the x algorithm is given by 



D. Non-local Actions and Noisy Hybrid Algorithms 



There is considerable interest in having the number of fermion flavours unequal to that described by the fermion 
kernel (e.g., less than four flavours for staggered fermions). For such theories the required probability distribution is 
given by P{U) oc e^^'^^'^^ detAl ([/)", where the number of multiplets n determines the number of fermion flavours 
(e.g., n = for Wilson fermions and n — \Nt for staggered fermions). For n ^ Z, the conventional pseudofermion 
approach fails because neither a non-integer power of the Dirac operator nor its derivative can be evaluated directly, 
which would be required to calculate the force. 



1. Ro Algorithm 



An alternative to the pseudofermion approach is to rewrite the determinant in trace log form where the effective 
fermion action is Sp = — ntrlnA^. The pseudofermionic force is replaced by a noisy estimator for the trace, since 
computing the trace exactly is prohibitively expensive. This force is written as 



dU 



-n tr 



-ntr 



= —71 tr 



d\nM 



dU 
(M'^M) 







= —ntr 


r du\ 



dU 



dU 



oU 



where 77 is a complex noise vector sampled from a Gaussian heatbath of unit variance. Defining an auxilliary field 
X = M'^-q the force becomes 



dU 



(11) 



With this formulation we can use the noisy Hybrid algorithm of f|V]to generate gauge field configurations corresponding 



to any number of flavours: the noisy estimator for the force S' being deflned by 
algorithm of Gottlieb et al. [10] and has leading order error 0{t5t). 



dU 



-n(E') 



This is the Rn 



Similar to the $ algorithm except that it grows with t because there is no shadow Hamiltonian. 
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E. Reversibility, Area Preservation, and the R Algorithm 

As can be seen from equations ([9]), p^ . and pT|) . the <&, and i?o algorithms have similar "pseudofermion" 
force terms despite their different derivations. Indeed, we introduced the x algorithm to emphasise this similarity: in 
both the X ^-nd. i?o algorithms the "pseudofermion" force is computed from Gaussian noise 77 for each MD step. The 

difference is that in the former the pseudofermion field x — \U{t)\ri is calculated from a heatbath at the beginning 



of each MD step, whereas in the latter the auxilliary field x — M'^ \U{t + ^5T)j 77 is calculated at the midpoint of each 

MD step (that is, at the same time as the force itself is evaluated for the integrator). 

The X algorithm has 0{t5t^) errors for n = 1 multiplets, whereas the i?o algorithm has 0{t5t) errors. However, 
for ri = multiplets (i.e., no fermions) the algorithms are identical and have errors of 0{5t^). It seems reasonable to 
expect that the leading error has a linear dependence both on the time within the MD step at which the pseudofermions 
are generated from their heatbath and on the number of multiplets, so if we evaluate the pseudofermion field at time 
t = i(l — n)(5T through the MD step for < n < 1 we should obtain an 0{t5t'^) algorithm. This is the R algorithm 
of Gottlieb et al. p^ . For two flavours of staggered fermions, this means evaluating the pseudofermion fleld a quarter 
way through each MD update. Note that this algorithm is neither reversible nor area-preserving. 

To prove that the R algorithm does indeed have 0{t5t'') leading order error we again look at the fixed point 
(equilibrium) distribution. The condition is that given in equation ([4]), and in this case we have neither area- 
preservation nor reversibility. We consider a single step of the R algorithm, where the auxilliary field x is computed 
at a time t = ^(1 — a)5T with a a free parameter. Expanding equation ([4]) to leading non- vanishing order we obtain 

{{5 + mS)^^^ ^-{{5 + 5)H + trln U.)^^^ + ■■■. 
We can compute the leading contributions to this quantity as follows; the change in energy over a trajectory is 



{5H)p^^ = 2n{n + 2p^a) tr 



dU 



dU 



where the 0{5t) term vanishes upon noise averaging. Taylor expanding the Jacobian for each update step gives 

(tr In J7,)p,,, = natr 



M ^—M 



dU 



dU 



X 2 

or . 



Finally, the leading order contribution to the quantity 6U = U ^{St) — FoUoF that measures the lack of reversibility 
of the integrator is 



{SH)p^,j — 2p^na tr 



M -— -A^ 



dU 



dU 



5t^ 



We thus find that 



where 



A = ^n{n - a)(l - ]?) tr 



M -— -A^ 



dU 



dU 



If we choose a — n the leading term cancels, and thus the leading error is 0(r(5T^) for an entire trajectory of t/St 
steps. Therefore, as claimed, the R algorithm has errors of 0{tSt'^), and thus so does the x algorithm since it 
corresponds to the special case of the R algorithm with n = 1. 



F. Exact Noisy Algorithms 

It is interesting to consider whether the noisy algorithms described above can be made exact (in the sense of having 
no integrator step-size errors) , and if so how. 

With the addition of an accept/reject step after the MD step, the $ algorithm becomes the exact Hybrid Monte 
Carlo (HMC) algorithm 0. 
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The noisy Hybrid algorithm can be made exact by including a noisy acceptance step [H, [s^, [H, [1^ after each MD 
integration step. Note that a trajectory is defined as being the MD evolution between momentum refreshments, and 
can consist of any number of MDMC steps, where an MDMC step is an MD step followed by a (noisy) acceptance 
test. Using just one acceptance test at the end of the molecular dynamics trajectory is not valid, since reversibility 
is violated because of the noise. Unfortunately this exact version of the algorithm suffers from a significantly longer 
autocorrelation time: this is because the momentum must be flipped if a rejection occurs at any of the t/6t noisy 
acceptance tests required per trajectory, just as for the second-order Langevin (Kramer's) algorithm 0, So, |4l|. 

The X algorithm can be made exact by the addition of a Metropolis acceptance step after each MD step, and 
including a momentum flip after every rejected MD update: however, the resulting algorithm suffers from the same 
problems as the exact noisy Hybrid algorithm discussed before. 

The Rq algorithm can be made exact through the addition of a noisy acceptance test, however, since the algorithm 
has scaling 0{tSt), this would require very small step sizes. Unfortunately the R algorithm cannot be made exact by 
adding a (noisy) Metropolis acceptance step, because the lack of area-preservation and reversibility preclude detail 
balance being satisfied. 

Clearly, the use of exact algorithms is preferable to that of the inexact algorithms described in this work. For non- 
local actions (e.g., a non- integral number of fermion multiplets) the Rational Hybrid Monte Carlo algorithm .lli[l2| is a 
good candidate: in this algorithm the fractional power of the fermion kernel that appears in the pseudofermion bilinear 
is replaced by a rational approximation that can be directly evaluated and differentiated. Hence a pseudofermion 
formulation can be used, and a Metropolis acceptance test can be added to render the algorithm exact. 

VI. CONCLUSIONS 

In this paper we have described how the use of a symplectic integrator for the numerical integration of Hamilton's 
equations of motion not only preserves the reversibility and area-preservating properties of their exact solution, but 
also conserves a "shadow" Hamiltonian close to the original one. We have shown how this conserved Hamiltonian may 
be computed as a power series in the integration step size St using the BCH formula for the Lie algebra of Poisson 
brackets. 

We then considered Markov processes of the Hybrid type, in which molecular dynamics and momentum refreshment 
Markov steps are alternated. Since the fixed point of the molecular dynamics step does not coincide with that of the 
momentum heatbath neither can be the fixed point of the full Markov process. We have derived a general condition 
(equation ^) for this fixed point distribution, which simplifies to equation ^ for the case of reversible and area 
preserving MD integrators, and shown how properties of this distribution can be found by expanding these conditions 
in powers of the integration step size St. For the case of the Langevin algorithm, it was shown how the equilibrium 
distribution can be found explicitly to any order in 6t, and why this is only an asymptotic expansion. 

Finally we considered those algorithms which use a noisy estimate of the force. Here the leading order behaviour 
of these algorithms was found for the noisy Langevin and Hybrid algorithms. It was shown that in general noisy 
Hybrid algorithms have leading order error 0(t(5t) regardless of the order of the numerical integrator, however, there 
are special cases where we can cancel this leading term through a judicious choice of when we evaluate the noise, 
i.e., the R algorithm. We also considered how to render these algorithms exact through the addition of a Metropolis 
acceptance test. 
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In fact, it is essentially the same algorithm. 
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APPENDIX A: HAMILTONIAN VECTOR FIELDS AND POISSON BRACKETS 

We shall denote by A'^ the set of antisymmetric multilinear fc-forms that act on fc-tuples of vectors in the tangent 
bundle TA4 over a manifold A4. 

Definition 1 The exterior derivative : A'^ — > A'^+^ is the unique linear transformation satisfying 

1. df{v) = v{f) for any 0-form f G and vector field v G TA4, 

2. d^ = 0, and 

3. d{a A /3) = (da) A /3 + {-l)'^^^°'a A dp (the anti-Leibniz rule). 
Lemma 1 For any O E A'^ 

dn{vo, ...,vk) = ^{-yvin{vo, ...,Vi,...,vk) 

i 

+ ^{-y'^^fl{[Vi,Vj],Vo, ...,Vi,...,Vj,..., Vk), 
i<j 

where Vj indicates that the term vj is omitted, and [a, b] G TA4 is the commutator of two vector fields a,b G TA4. 
Proof. We may express f2 in a local coordinate patch as O = ^fli^^^,,,^ij,^dq'^^ A ... A dq^'', and thus 

dQ{vo, ...,Vk) = ^ ^"g^;;^" ^?^" A dq^^ A ... A dg« (^;o, . . . , Vk) 
= ^(_)i ^^W....,A.,-.W ^^o^Mi . . .^M. . . .^/-^ 



i=0 



Vn ■ ■ ■ vy ■ ■ -v 



J 



^1 

k 

dq^' 

o<i<j<k L ^ ^ . 

= ^{-yViCl{vo, ...,Vi,...,Vk) + '^{-y~^^fl{[Vi, Vj],Vo,...,Vi,...,Vj,...,Vk). 

i i<j 

The result is independent of the coordinate system used for this verification. ■ 
For any 1-form and 2-form uj this identity is 

d0{a,b) = a0{b)-be{a)-0{[a,b]y, 
dii){a,b,c) = aw(6, c) + 5cj(c, a) + ccj(a, 6) (Al) 
-oj{[a,b],c} -uj{[b,c],a) - uj{[c,a],b). 

Definition 2 The cotangent bundle T*Ai has a symplectic structure if there is a non-singular closed fundamental 
2-form uj. 

Definition 3 For each 0-form F on T*M there is a corresponding Hamiltonian vector field F defined by dF = ipuj 
where i is the interior produce; equivalently we may write this as dF{x) = u{F,x) Mx G TA4. 

Definition 4 The Poisson bracket of two 0-forms is 

{A, B} = -Lo{A, B) A,BgA°. 

Lemma 2 The action of a Hamltonian vector field A on a 0-form F is given by AF = {A, F}. 

Proof AF = dF{A) = ipOj{A) = uj{F, A) = {A, F}. ■ 
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Lemma 3 The space of 0-forms on T*A4 together with the Poisson bracket form a Lie algebra; that is the Poisson 
bracket satisfies 

• {A, ^} = and 

. {A {B, C}} + {B, {C, A}} + {C, {A, B}} ^ 
(Jacobi identity). 

Proof. The antisymmetry is obvious. To establish the Jacobi identity consider duj(A, B,C) for three arbitrary 
Hamiltonian vector fields A, B and C. Recalling that the fundamental 2-form is closed, duj = 0, and using equa- 
tion (jAip we have 

duj{A,B,C) = Auj{B,C)+I3uj{C,A) + Cu{A,B) 
+Lui[A, B],C)+uj{[B,C],A) 
+u;{[C,A],B) 

= 0. (A2) 

Now, Aijj{B, C) = —A{B, C} by the definition of the Poisson bracket, so Auj{B, C) = —{A, {B, C}} by application of 
lemma [5] 
Similarly, 

u;{[A, B],C) = -u;{C, [A, B]) = -i^ojUA, B]) 
= -dC{[A,B]) = -[A,I3]C 
= -{AB-BA)C = -A{B,C} + B{A,C} 
= -{A,{B,C}} + {B,{A,C}}. 

The Jacobi identity follows upon subsitituting these relations into equation IA2I ■ 

Theorem 1 Let A and B be two Hamiltonian vector fields, then their commutator is also a Hamiltonian vector field 
and furthermore [A, B] = {A,B}. 

Proof. Consider the action of the commutator on an arbitrary F £ A", [A, B]F = {AB — BA)F = {A, {B, F}} — 
{B, {A, F}} upon applying lemmaH Using the Jacobi identity we obtain [A, B]F = -{F, {A, B}} = {{A, B}, F} = 
{A, B}F. As this must hold VF the result follows. ■ 

This argument may be carried out more explicitly in a local coordinate patch. Locally u may always be written as 
dqi A dpi (Darboux theorem) so 



uj{A,x) = dA{x) 



dqi A dpi^ ^a- 



id 9 u d ,u d 

1 --fa-' r X \- X 

dq^ dp^ ' dq^ dp^ 

^~^ ~. , dA dA^^ 

ax —ax = TT-x + -TT-^x ; 

oq^ op^ 



and as this must hold for arbitrary x we may identify 

^ OA ^ dA 
op^ oq^ 

and thus 

'dA d dA d 



dq^ dp^ 
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The commutator of A and B is 



[A,B] 



E 

i 

= E 



+ 



dA a 

dA d'^B 
dp^ dq^dpi 

dA d^B 



dA d 
dq^ dp^ 



dB d'^A 
dp^ dq^dpi 

dB d'^A 



dB d 
dpi dqi 

dA d'^B 
dq^ dp^dpi 

dA d'^B 



dB d \ 

dqiW) 



+ 



dB__8^A_ 
dq^ dp^dpi 

dB d^A 



dp^ dq^dqi dp^ dq^dqi dq^ dp^dqi dq^ dp^dqi 



d 

dqi 

d 

dpi 



= E 



-y 

dpi ^ 



( d{A,B} d _ d{A,B} d 



dA dB 


dA dB\ 




d 


fdA dB 


dA dB\ 




Qpi Qqi 


dq^ dp^ J 


dqi 


dqi 


ydp^ dq^ 


dq^ dp^ J 


dp) 



E \^ dpj dqi 



dqi dpi 



{A,B}, 



where the Poisson bracket is 



dA dB 
dq'- dp^ 




dA dB 

Qpi Qqi 



dA d 
dpi dqi 



dA d dB d 
dqi dpi ' dp^ dq^' 



dB d 

Qqk Q.pk 



APPENDIX B: BASIC PROPERTIES OF LIE ALGEBRAS 

Definition 5 Let K be a commutative ring with unit. A Lie Algebra over K is a K-module^^ L together with a 
K-bilinear mapping L x L ^ L : {x,y) i-^ [x, y] called a Lie bracket which satisifies 



[x, a;] = 

[a;, [y, z]] + [y, [z, x]] + [z, [x, y]] = (Jacobi identity) 



(Bl) 



for all x,y,z G L. 

Note that this implies that the Lie bracket is antisymmetric, 

[x, y] + [y, x]=0 



Vx, t/ e L 



because [x+y, x+y] — [x , x] + [x , y] + [y , x] + [y , y] = [x, y] + [y, x] . The converse is also true unless K has characteristic 2. 
If A is an associative algebra over K, then it has a natural Lie algebra structure^^ given by [x, y] = xy — yx. 

Lemma 4 For any given Lie algebra L there is a unique associative algebra Ao called the enveloping algebra of L with 
a Lie algebra homomorphism 4>q : L ^ Ao defined by ipQ : [x, y] i-^ xy — yx which has the universal property that for 
all A for which there is a Lie algebra homomorphism (j) : L^ A there is a unique algebra homomorphism f : Ao ^ A 
such that <p = f o <pQ. 

In other words the enveloping algebra "contains" all associative algebras which have L as their natural Lie algebra. 

Proof. Let T be the tensor algebra^^ over L and I be the ideal of T generated by the elements x®y — y®x — 
[x, y] (Va;, y e L); then Ao=T /I. ■ 



Recall that if X is a field then a JsT-module is a linear space. 
Note that the Jacobi identity holds automatically. 

T is the universal algebra over the module L. We shall denote multiplication in T by ®. 
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1. 



Poincare—Birkhoff— Witt Theorem 



Not only is a Lie algebra contained in its enveloping algebra, but in fact the tensor algebra has a direct sum 
decomposition T = 7^ ®X, where 7^ is the subalgebra a T consisting of all symmetric tensors, so — Ts- This is 
the content of the following 

Theorem 2 (Poincare— Birkhoff— Witt) Let L he a Lie algebra over K which is a free K-module with a totally 
ordered basis {xi), and let Ao be its enveloping algebra. Then Ao is a free K-module with the set of ordered products 
(t)o{xi^) . . .(/)o(xi,J (io < • • • < in) as a basis. 

Proof. Recall from the proof of Lemma 2] that — T/T. If we define Tg to be the submodule of T spanned by 
the ordered products Xi^ ■ ■ ■ Xi^ (ii < • • • < then what we must show is that T is the direct sum of modules 
Tg and T, that is T = 7^ © X and Tg OX — 0. We shall do this by showing that each element of T has a unique 
decomposition into an element of 7^ and an element of the ideal X. Clearly, for elements of Tg or X we must have 



a 1^ {x 1^ y — y 1^ X ~ [x,y]) 1^ P 
= ® ig) {x ig) y - y X - [x, y]) <S) /3|(a,/3 G T). 

For any other basis element of T where the factors are not in increasing order we have 

a g) y g) X g) P 

= |a (g) X (g) y ® /3 — a g) [a;, y] (g) /?! 

Q)lag){yg)x-xg)y+[x,y])g)Py 



where x < y. All that remains to show is that this decomposition is well-defined and does not depend upon the order 
in which we apply the preceding identity. The only non-trivial case occurs when there are three adjacent out-of-order 
factors, for which we have both 



agzgygxgp — lag)zig)xg)yig)l3 — ag)zg)[x,y]ig)(3>Q) <aig)zig){yig)x — xg)y+[x,y])ig)f3 



and 

a(g> z (g> y g X g P = lagygzgxgp — a'g[y,z](gx(gp>(B \a (g {z (g y — y g z + [y, z]) (g> x (g> P 






+a g) {z g) y — y g) z + [y, z]) g) X g) P 
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These two values differ by 

|a (g) a: (g) [?/, z] (g) /3 + a (g) [a;, z] (g) y (g) /? + a (g) z (g) [a:, y] (g) /? 

—a (g [a;, (g z (g /3 — a (g ?/ (g [x, z] (g /? — a (g [y, z] (g a: g) /?! 
© |a g) a; g) [y, z] g) /3 + a g) [a;, z] g) ?/ g) /? + a g) z g) [x, y] g) 
—a g) [x, g) z g) /? — a g) y g) [x, z] g) — a g) [y, z] g) a; g) 
= |a (g (a; (g [y, z] - [y, z] (g a;) (g ^9 + a (g ([a;, z] (g y - y (g [a;, z]) (g /? + a (g (z (g [a;, y] - [a;, y] (g z) (g (3^ 

© |a ® [a;, z] g) y g) /3 + a g) z g) [a;, y] g) /? — a g) [a::, y] g) z g) — a ® y g) [a;, z] g) /3 — a ® [y, z] g) a; g) /3 
= {a g) ([a;, [y, z]] + [y, [z, x]] + [z, [x, y]]) g) /?} © |a g) ([x, [y, z]] + [y, [z, x]] + [z, [x, y]]) g) /jj = 
where we have used the Jacobi identity. ■ 



2. Free Lie Algebras 



We are interested in the properties of the Lie algebra generated by some set of operators, but a priori we know 
nothing about the nature of the Lie brackets of these generators. We therefore wish to work in the context of the 
most general Lie algebra which can be constructed from these generators, for which all Lie brackets are assumed 
distinguishable unless they are related by the defining relations (jBl[) : any further relations between Lie brackets may 
be applied post facto. More formally this means that wc wish to carry out our calculations in the free Lie algebra over 
our set of generators. 

Let Lq be a Lie Algebra over K, and A a set with a mapping i : A ^ Lq. Lq is called free on A if for any Lie 
algebra L and any mapping f : A ^ h there is a unique Lie algebra homomorphism / : Lq ^ L such that f o i ~ f . 
This is a universal definition, as shown by 

Theorem 3 For every set A there is a unique free Lie algebra L{A) on A, L{A) is a naturally graded K -module, i is 
an injection, the component of L(A) of degree 1 is the free submodule generated by i{A), and L[A) is generated as a 
Lie algebra by A. 



3. Hall Bases 



Central to calculations in free Lie algebras is the question of how many independent basis elements are there of a 
given degree, and how to reduce an arbitrary expression to canonical form in terms of such a basis. 

Definition 6 (Hall Trees [42], |43| . l44l |45| |) Given a set A consider the set of all binary trees with leaf nodes labelled 
by elements of A: this set is called the free Magma M{A). We shall denote the tree h whose left subtree is h' and 
whose right subtree is h" by h — [h' , h"). A Hall set H is a totally ordered subset of M{A) containing A which satisfies 

h<h" yh= {h',h") e H - A 

h = {h',h") e H 
h', h" e H and h' < h" and 
either h' Cz A or h' — (x, y) and y > h" . 

There is map / : M{A) L{A) defined by /(a) ^ a if a e A, and f : {h' , h") ^ [f{h'), f{h")]; the result of applying 
this map to a Hall tree gives a Hall word, and the Hall words corresponding to any Hall set form a Hall basis for L(^) 
(as a module). 
We shall use the following Hall basis for our calculations: x G iff 

X e A 

or X = [y, z] y,z e H, y < z 

or x=[y, [z,u]] y,z,ueH,z<y 

and X < y if degx < degy. 
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Any expression built out of Lie brackets may be reduced to canonical form by the applying the following transforma- 
tions: 

[sx + ty,z] 1-^ s[x , y] + t[y , z] where s,t G K 

[z,y] 1-^ -[y,z] iiy < z 

[y, [z, u]] i-> [y, z]] + [z, [y, u]] ii y < z < u, 

with the elements of A themselves ordered lexicographically (i.e., alphabetically). 
For generating set {A, B} this leads to the following basis: 

{a b}, {[a,b]}, 

{[A,[A,B]l [B,[A,B]]}, 

{[A,[A,[A,B]]], [B,[A,[A,B]]], [B,[B,[A,B]]]], 

{[[A,Bl[B,[A,B]]], [B,[B,[B,[A,B]]]l 

[A,[A,[A,[A,Bm, [B,[A,[A,[A,Bm, 

[[A,Bl[A,[A,B]]], [B,[B,[A,[A,B]]]]}, 



4. Dimension of Hall Bases 

The dimension of the Hall basis of degree iV on a set of cardinality q is given by Witt's formula [i^ ajv = 
Jr Sd|jv t^{d)q^/'^, where is the Mobius function. 

By the Poincare-Birkhoff~Witt theorem the set of ordered (symmetric) monomials on the independent commutators 
is a basis for the universal enveloping algebra of the free Lie algebra. A basis for words of length N is therefore 
provided by symmetric products of n^, words of length k chosen from the generators of the free Lie algebra, where 
Efc>i ^"-fc = ^- There are exactly (—)"'' (^°''') — C'''''^^''"^) ways of choosing these n^. words symetrically, and this 
is the coefficient of in the series expansion of (1 — x^)^"-'' = X]nfc>o {^nk)i~^'^)"'° ■ total number of words 

of length N is thus the coefficient of x^ in the generating function g = nfe>i(-'^ ^ x'')^"'''. On the other hand, the 
universal enveloping algebra is just the free algebra on q symbols, so there are q^ independent basis elements for 
words of length TV, and thus is determined from the equation g = J2T=ol'^^'^ = (1 ~ Qx)~^. Witt's solution is 
obtained by taking the logarithm of this equation, — J2k>i ln(l — x^) = — ln(l — qx), and equating the coefficients 
of x^ , ^d\N o-d/{N/d) — q^ /N. Using the Mobius inversion formula we obtain aN ~ J2d\N t^{d)q^^'^- 

The number of independent commutators on q letters is therefore 

^9(9-1), 
\q{q-l){q + l), 
\q\q-l){q + l), 
\q{q-l){q+l){q^ + l), 
i(7((7-l)(g+l)(g3 + g-l), 

More specifically we have o!f^ = 2, 03^^ — 1, ag^^ — 2, a-^'' — 3, a'^^ — 6, Og^' — 9, a'^^ = 18, Og^' = 30, a'^^ = 56, 
a'^iQ — 99, a''ii — 186, 0^3^ = 335, and so on. 



4^^ 
4^^ 
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5. Baker— Campbell— Hausdorff Formula 



Theorem 4 ([42]) Let M be an associative algebra over a commutative ring K Q and let d be a derivation on M. 
For any power series f{X) ~ X]ri>o '^n^^ we have 



k>l 



where ad a : 6 i— > [a, 6] . 

Proof. Equation (IB2p is linear in /, so it suffices to consider f{X)= X", for which we have 



k=l 



We shall prove this by induction on n: for n = it is trivially true, and 



(B2) 



k=l 



= d{X X") ^dXX'^ + X dX" ^dXX" + Xj2(y) (ad X)''-\dX)X' 
Using the identity Xu = [X, u] + uX = (ad X)u + uX we obtain 

= dXX" + ^ {&dXf{dX)X''-'' + Y^ ( ")(adX)'=-l(dX)X"+l- 



-/c 



k=l 



k=l 



n+1 



fc=2 



X" + ^ f J (ad X)''-\dX)X''+^-''- + (ad 



^dxx'' + J2 



k=2 



n 
k-1 



k=l 

\k — lr^iv\ vn+l-k 



{a.dXf-'idX)X''+'-'' + (adX)"(dX) +ndXX'^ 



n+l 



fc=l 



{adXy-\dX)X''' 



Corollary 1 



de^ = g{sidX){dX)e^ 



where 



1 _ 1 



/c>l 



k\ X 



Proof. By Theorem 2] we have 



de^ = E ^(adX)'=-i(dX)e^ = 5(adX)(dX)e^ 



fc>i 



Definition 7 We define the Hausdorff series as 

H = ^Cn{A,B) wh 



ere e" = e^e^. 



n>l 



and the c„ are homogeneous of degree n in the generators A and B. 
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Theorem 5 Using the derivations^^ (Ia and ds defined by dAA = A, dAB = 0, dsA = 0, and dsB = B , and setting 
d = dA + ds, we have 

dH = i^-^ coth ^ j (A + B)- i^-^j {A - B). 

For any derivation D we have Dl = D{1 ■ 1) = Dl • 1 + 1 • Dl = 2 so D\ = 0. It is also trivial to verify inductively 
that dxX"" = nX" and d^e-^ = Xe^ . 
Proof. Using Corollary [1] we obtain 

dAc" = dAie'^e^) = dAc"^ + dAC^ 



dse'" 



whence 



and 



Ae" = g{a.dH){dAH)e" , 

dB{e-^e~^) = dee^^ e"^ + e'^ dge-^ 

-Be-" = g(ad(-i/)) e"^; 



A = giadH){dAH), 
B ^ g{-&dH){dBH) 



dAH = [g[adH)]-^A, 
dsH = [g{-iidH)]-^B. 

We now observe that 1/5(0;) may be decomposed into the sum of an even and an odd function 

1 X / X \ 

— 7 — — coth — =F 1 , 

g{±x) 2 V 2 ^ / ' 

and the desired result follows immediately since d ~ dA + ds ■ ■ 



Definition 8 The Bernoulli numbers B„ are defined by 

1 X 



g{x) 

Taking the symmetric part of this expression we see that 



X , X \ - 

— coth — = > 
9 9 ^ 



B^x'- 



(2m)! 



Theorem 6 (|47l. l48l . l49l . l50l . l5lL l52l |) The terms in the Hausdorff series Cn{A,B) are given by the recursion 
relations 



Cn+l 



^A-h[Cn. ^ - 5] + E T^^i E [^^- [■ • ■ ' [^^-' ^ + B]..]]\. 

fclH hfe2m = i 

Proof. Using the properties that the c„ are homogeneous in A and B and that ad 1 = we have 

dH = ^dcn ^^ncn, 



n>0 



n>l 



adH — E] s-d c/c ; 



fc>i 



15 If = R these definitions are equivalent to d^/{A, B) = ^-^^^f-^^ I , ds/(A, _B) = ^f'-'^f^^ I and B) = I 
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dH = I coth^]{A + B)- (^ \ (A-B) 



and so 



hence by Theorem [5] 

riff = , ^„„.. , , ^ , 

V 2 2 y ' V 2 

m>0 ^ ^ ^ ^ 

m>0 ^ ' k>l ^ ^ 

L"/2J „ 

^(n+l)c„+i = E E 7^ E adcfc,...adc,,^(A + i3)-^i(adc„)(A-i3) 

n>0 n>l m=0 ^ '' fci , • - • .fc2„i > i n>l 

fcl + --- + fe2m = " 

therefore, equating terms of equal gradation, 

(n+l)c„+i ,^ "', adcfci . . . adcfc2,„(A + S) - i(adc„)(yl - B). 

■ 

The first few terms in the HausdorfF series are 

In(e^e^) = {A + S} + S] + ^{ [A B]] - [S, [A S]]} 
[AS]]] 

+^{ -4 [S,[A[A[Ai3]]]] -6 P,i3],[A,[A,i?]]] 

+4 [B,[i?,[A,[A,i3]]]] -2 [[A,S],[B,[A,S]]] (B3) 
- \AAAAAAA,B\\\\ + [S,[i3,[S,[AS]]]] } + ••• 

From this we easily obtain the formula for a symmetric product 

ln(e-^/2e^e^/2) = {A + S} - Jj{2[i3, [A, S]] + [A [A S]]} 

+ 5=L_| 32 [S,[i?,[A,[Ai3]]]] -16 [[A,i3],[B,[A,i3]]] 

+28 [B,[A,[A,[Ai3]]]] +12 [[AB],[A[A,i3]]] (B4) 
+8 \B,\B,\B,\A,B\\\\ +7 [A, [A, [A [A, B]]]] } + ••• 

Since e^e^ = gA+B+5(A,i3) ^^^^^ g-^e"^ = g-s-A+5(-B,-A) ^j^^^^ = -(5(-B,-A), so under inter- 

change of A and B all terms of even grading in change sign, whereas those of odd grading do not. 



APPENDIX C: HIGHER-ORDER SYMMETRIC SYMPLECTIC INTEGRATORS 

It was observed by Campostrini [l^,!!^ that one can construct higher-order integrators by the following method p^ : 
since l!<i(5r) = e^''^ + Ba^St'^ + 0((5r^) we observe that the "wiggle" 

C/o(e)[/o(-ae)!7o(e) = e^^^-.)^ ^ ^^^3 - a3)e3 + 0(,5t^) 

can be adjusted to give an integration scheme correct to O(^r^) by choosing a — \f2. The step size may be kept fixed 
by taking e = 8t j {2 — a). 

Naturally, this wiggle may itself be iterated to give integration schemes of arbitrarily high order, all of which are 
still reversible and area-preserving. We use the recursive definition 

Unie„)Un{-a„en)Unien) 
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and choose (t„ = ^"+\/2 and = (5t/(2 — (t„). 

It was noted by Yoshida [l6| that while the lowest-order Campostrini scheme is optimal, in the sense of requiring 
the fewest integration steps, the second order Campostrini wiggle is not. The lowest order Campostrini scheme 
corresponds to the operator 

g€Q/2gePgeQ/2g-<TeQ/2g-<T£Pg-aeQ/2g€Q/2g£Pg€Q/2 

^ geQ/2g£Pg(l-<T)eQ/2g-rTePg(l-<T)£Q/2gePg£Q/2 

which uses precisely seven steps. The second order wiggle uses 19 steps, while the 15 step operator 

Uo{u St) Uo{v 6t) Uo{w 5t) Uo{x 6t) 
Uo{w5t)Ug{v 6t)Uo{uSt) 

_ gU StQ/2^u StP^(ii+v) StQ/2^v StP ^{v+w) StQ/2 
gU) StP ^(w+x) StQ/2^x StP ^{w+x) StQ/2 
gU) StP ^{v+w) 5rQ/2gD StP ^{u+v) 5tQ/2 
^uStP^u StQ/2 

has errors of 0((5r^) if the ideal defined by 

+ 2w^ + 2v^ + 2u^ = 0, 
+ 2w^ + 2v^ + 2u^ = 0, 
x + 2w + 2v + 2u = 1, 
ixwv — 2xw V — 2x vu -\- Axwu — Aw vu 
— 2x — 2xw u — 2x wv + 8wvu — 2xv u 
+4:xvu^ ~ x^w + x^w^ - Awv^u - 2w^v'^ 
-2v^u^ - 2w^v? + 2wu'^ + 2vu^ + Av^u^ 
-Aw^u + iw'^u^ ~ Av^u + 2wv'^ - x^u^ 
- x'^u + x^%i^ - Aw'^v + 4w^w^ 
-\-x V — X V + XV — X V — X w + xw — 

is not empty. This is indeed the case, as a Grobner basis computation shows that w, w, and x may be expressed as 
polynomials in u, which is a root of the following irreducible polynomial 



5632424294400000000^39 - 922003368192000000001*^8 
+7106323614105600000001*37 - 3437629764814080000000^36 
+11745037928943360000000^35 - 30260229452421120000000^34 
+61344468339328512000000^33 - 100894346480650176000000^3- 
+137871545973425856000000^31 _ 159597428255349696000000m' 
+159057766014056179200000^29 - 138323253491741289600000m' 
+10609941741061 1328000000^27 - 723618101 16050054400000^^* 
+44125123305044761920000u25 _ 24138408506765309280000^^^ 
+11867110408796028480000^23 - 52470119653215278400001*22 
+20868001525237579200001*21 - 746466059135744064000w20 
+2401102666276079040001*1'^ - 694318771425474720001*1^ 
+18041618760883056000u" _ 4210061488653312000^^6 
+8814266341001568001*15 - 1653355743058944001*1'^ 
+277318847797704001*13 - 4148250096765600^12 
+551410054740000^11 - 64829840769360^1° 
+6700635295200u^ - 6040164600001*^ 
+46995575760^7 - 3112757280^^ 
+172255032^5 _ 7756920i*'i 
+273360^3 - 70801*2 
+120u- 1. 
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This polynomial has three real roots, corresponding to the numerical solutions found by Yoshida: 



11 = 


0.784513610477557263819497633866, 


V 


= 0.235573213359358133684793198233, 


w = 


-1.177679984178871006946415596562, 


X 


= 1.315186320683911218884249687935; 


u = 


1.439848167976783090930499281479, 


V 


= 0.004260681870792016799146793392, 


w = 


-2.132285222001451515523597811357, 


X 


= 2.376352744307752823717294656042; 


u = 


1.447782562399297932897896663298, 


V 


= -2.144035316305389310213622103136, 


w = 


0.001528862284249274922787428878, 


X 


= 2.389447783243684212186399178641. 
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